Comment on "Bose-Einstein condensation with magnetic dipole-dipole forces" 
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The ground state solutions of a dilute Bose condensate with contact and magnetic dipole- dipole 

interactions a re examined. By lowering the value of the scattering length, Goral et al. [ cond- 

mat/9907308| and Phys. Rev. A 61, 051601 (2000) ] numerically predict a region of unstable 
solutions, accompanied by a neighborhood where the ground-state wave functions have internal 
structure. On the contrary, we find that the dipolar condensate has an intuitively-located stability 
region, and ground-state solutions near the instability threshold which are absent any unusual 
structure. 

PACS number(s): 03.75.Fi, 32.80.Pj, 03.65.-w 

In a recent paper Q, Goral et al. studied a non-ideal 
Bose-Einstein condensate in the presence of magnetic 
dipole-dipole forces. For a magnetic moment equal to 
that for chromium and a typical value for the scatter- 
ing length, all solutions were found stable and differed 
only in size from condensates lacking long-range interac- 
tions. Upon reducing the scattering length below a criti- 
cal value they observed the dipolar condensate to become 
unstable, and for scattering length values just above the 
instability a "structured" condensate was found (i.e., the 
density had several peaks instead of the ordinary single- 
peaked behavior) . The purpose of this Comment is to re- 
port a disagreement in both the location of the instability 
threshold, and the existence of the structured solutions. 

The results of Ref. § were achieved by numerically 
solving the Gross-Pitaevskii equation for atoms in a 
cylindrical harmonic trap: 



where a is the angle between k and fi. 

Now, inspecting the formula (^), one thing seems in- 
tuitively clear: the first term will introduce an effective 
shift in the scattering length according to 



U = 



(4) 



and one would expect instability (roughly) when this be- 
comes negative. A first approximation to the instability 
threshold is then 



ac ~ 12nh 2 • 

Using the chromium parameters fi 
get the estimate (aNa = 2.75 nm) 
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Here '3/ is the condensate wave function, a is the s-wave 
scattering length, N is the number of atoms, and m is 
the mass of the atom. The long-range potential is due to 
the magnetic dipole-dipole interaction and is given by 
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where u = (r — r')/|r — r'| and is the magnetic per- 
meability of the vacuum. All the magnetic moments are 
assumed to point in the same direction (^-direction), i.e., 
fii = fi2 = fJ-z. As explained in Ref. |Q, the Gross- 
Pitaevskii equation (Q) can be solved easily enough when 
one notices that the nonlocal term is a convolution, and 
is therefore local in fc-space. This approach requires a 
Fourier transform of the dipole-dipole interaction, given 
in the limit of small atomic overlap distance as 



This simple estimate is about twice as big as the values 
given by Goral et al. in their Fig. 2, and it also indi- 
cates a critical scattering length that is independent of 
particle number. To a certain extent, the kinetic energy 
(and possibly the anisotropic part of the dipole-dipole 
interaction) can stabilize the condensate in a manner 
similar to that for condensates with negative scattering 
lengths, but it seems unlikely that these would cause such 
a dramatic shift in the critical scattering length. 

To determine more accurately when the condensate is 
expected to turn unstable, we perform a variational cal- 
culation by assuming a trial wave function for the con- 
densate density 



|*|2 = 
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■ exp (-(a; 2 + y 2 )/(J 2 ) exp (-z 2 /<r 2 z ) 
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The energy of the condensate is then E = Ek 
E NL + E DD , where 
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P[V(r)] = -i/i 0( u 2 (l - 3 cos 2 a), 
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Here 2Fi(a, (3; C; 2) is a hypergeometric function, and 
Eq. (||) has been used to calculate the anisotropic con- 
tribution (Edd) to the energy. Minimizing this en- 
ergy functional for a spherically symmetric trap with 
u) = (2tt)150Hz and N = 300 000, we find that the en- 
ergy becomes unbounded and the system unstable when 
(i/«Na ~ 0.94 a c . Additionally, the critical value is quite 
insensitive to the number of particles: for N = 10 7 it is 
a/aNa ~ 0.95 a c . These values are in disagreement with 
the results of Ref. §. 

We also used a trial wave function 



l*l 2 = 
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exp(-(z 2 + y 2 )/<7 2 ). (12) 



At the threshold of instability this trial wave function 
seems better justified than the Gaussian trial. Due to 
the cos 2 a factor in Eq. (||) , there is extra repulsion in 
the z-direction, so inverted parabola seems to be a natu- 
ral guess. In the xy-plane the interaction terms vanish so 
a Gaussian profile is expected. For the function (]lj) the 
kinetic energy diverges so it was left out from the opti- 
mization process. Thus this approach is valid only at the 
Thomas-Fermi limit, i.e., at large particle numbers. We 
were not able to evaluate all the integrals analytically, so 
a numerical integration was performed. With this wave 
function we estimate that for N = 300 000 the critical 
value is a/aNa ~ 0.96 a c , again quite close to the simple 
estimate given before. 

In Fig. [I] we show the instability threshold calculated 
with three different methods. It can be seen that each 
method gives about the same results at large parti- 
cle numbers. The missing kinetic-energy term for the 
parabolic-Gaussian trial results in a non-physical behav- 
ior at small particle numbers, but the Gaussian trial gives 
fairly satisfactory results for all particle numbers. Fig. |l| 
should be contrasted with Fig. 2 of Ref. |TJ . 

Finally, we address the existence of structured dipo- 
lar condensates by solving Eq. (jl]) numerically using the 
method outlined in Ref. Q| . The largest grid we used had 
a size 128 x 128 x 128, but the results obtained with this 
grid were not very different from those in a grid with a 
size 64 x 64 x 64. For N = 300 000 we estimate the critical 
value of the scattering length to be a/a^ a ss 0.92 a c . Near 
the instability threshold, the ground-state wave functions 



looked roughly parabolic along the z-axis and Gaussian 
in the xy-plane. Solutions in the neighborhood of the 
Ref. instability were consistently found to be unsta- 
ble, although transient structure was observed en route 
to collapse when the kinetic energy term was artificially 
removed. 

While it is true that a variational wave function will not 
deliver any more physics than is already present in the 
trial wave function, our analytical and numerical results 
strongly indicate that the numerical instability threshold 
of Ref. is dubious; furthermore, their structured so- 
lutions are not numerically reproducible. In conclusion, 
we therefore contend that, upon lowering the strength of 
the contact interaction, a non-ideal dipolar BEC reveals 
an instability threshold which is largely independent of 
N and given to a good approximation by the intuitive 
result (|6|), and the ground-state wave functions at and 
near this threshold have a simple, single-peak profile. 
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FIG. 1. Instability threshold as a function of the number 
of particles for a non-ideal dipolar condensate of chromium 
atoms in a spherically symmetric trap with uo = (27r)150Hz. 
The solid line is based on a Gaussian trial, and the dashed line 
is based on a parabolic-Gaussian trial. Individual data points 
are numerical solutions of the GP equation with error bars 
indicated. (The small amount of noise present in the figure is 
due to the Monte-Carlo method used for minimization.) 
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